Tutorial3

Analisis descriptivo de datos

#personalizar nuestra seccion

options(digits=3,papersize = "letter")  
# https://stat.ethz.ch/R-manual/R-devel/library/base/html/options.html


# Opciones globales en graficas
par(mar=c(5.1,5,4.1,2.1),font=3,family="sans",bg="yellow")   

1 Problema a estudiar

Una empresa quiere introducir un programa de estímulos para su personal. Debe hacerlo lo más justo posible para evitar disgustos y un ambiente pesado entre los empleados. La información que ha recabado es de una muestra aleatoria de 200 empleados.

• Antiguedad: antigüedad en años
• Hora extra: horas extra trabajadas en el año
• Sexo: 1-mujer, 2-hombre
• Cursos: número de cursos de capacitación en un año
• Incapacidad: días de incapacidad
• Aptitudes: índice de aptitudes y capacidades
• Escolaridad: 0-bachillerato, 1-licenciatura sin título, 2-licenciatura con título, 3-posgrado
• Salarios: salario anual en dólares
• Edad: edad en años del trabajador

1.1 Lectura y manejo de datos

Antes de realizar el análisis estadístico más adecuado es conveniente revisar:

1. Que la información de todas las variables sea completa
2. Que nuestra elección esté de acuerdo con el tipo de variable y su escala de medición para ¡no meter la pata!

setwd("C:/Users/Emman/Documents/CCH-N} FES-A/fessaa/FES/7mo semestre/Temas Selectos de Estadistica")
getwd()

## [1] "C:/Users/Emman/Documents/CCH-N} FES-A/fessaa/FES/7mo semestre/Temas Selectos de Estadistica"

datos1 <- read.table("./DATA/Empleados.txt",header=TRUE,row.names=1)

head(datos1)    # ver los primeros 6 renglones del conjunto de datos

##   antiguedad horasextra sexo cursos incapacidad aptitudes escolaridad salario
## 1         11        125    1      4           9     121.9           2   23065
## 2         24        225    2      2           2     114.2           1   27180
## 3         17        115    2      3           5     134.1           1   34875
## 4          9        117    1      1           1     114.0           1   23685
## 5         15         26    1      2           0     151.4           2   33550
## 6          6         43    1      4           3      96.7           1   22635
##   edad
## 1   44
## 2   50
## 3   48
## 4   53
## 5   62
## 6   45

tail(datos1)    # ver los ultimos 6 renglones del conjunto de datos

##     antiguedad horasextra sexo cursos incapacidad aptitudes escolaridad salario
## 195          5        173    2      2           8     122.1           2   29000
## 196          1         11    2      1           4      93.9           3   30200
## 197          7        129    2      3           6     104.6           1   34450
## 198          2        162    1      3           3     107.8           1   23550
## 199          2          5    2      1           1     101.7           1   22000
## 200          5         74    2      3           6     111.0           1   33000
##     edad
## 195   29
## 196   30
## 197   35
## 198   25
## 199   26
## 200   34

dim(datos1)

## [1] 200   9

names(datos1)

## [1] "antiguedad"  "horasextra"  "sexo"        "cursos"      "incapacidad"
## [6] "aptitudes"   "escolaridad" "salario"     "edad"

str(datos1)

## 'data.frame':    200 obs. of  9 variables:
##  $ antiguedad : int  11 24 17 9 15 6 4 2 17 17 ...
##  $ horasextra : int  125 225 115 117 26 43 124 71 166 158 ...
##  $ sexo       : int  1 2 2 1 1 1 2 2 2 1 ...
##  $ cursos     : int  4 2 3 1 2 4 2 1 2 3 ...
##  $ incapacidad: int  9 2 5 1 0 3 4 1 5 2 ...
##  $ aptitudes  : num  122 114 134 114 151 ...
##  $ escolaridad: int  2 1 1 1 2 1 2 1 1 1 ...
##  $ salario    : int  23065 27180 34875 23685 33550 22635 19575 20430 18955 25595 ...
##  $ edad       : int  44 50 48 53 62 45 26 28 33 40 ...

mode(datos1)

## [1] "list"

# Exploremos algunas variables

datos1$antiguedad

##   [1] 11 24 17  9 15  6  4  2 17 17 15 21  4 12 23 20 19 12  5 11 11  8 20  1  6
##  [26] 18 21  7 21 27 20 11 11  3 16  2 12 16  9 15  3 17 17 23  6  1  7  4 22 18
##  [51] 22 25 15 24 15  2 17  7  8 11 10 10  8 15  8 23 24 22  6  3 16  4 12 23 13
##  [76] 11  4  3  9 12  3 16 23  2 18  3  7 25  2 17 22  6 27 14 12 24 14 14 11  4
## [101]  3 12 12 19 12  4  5 23  9 24 21 19 14  3  3  8 15 18  7 11  5 18 12  2 26
## [126] 26 11 11  0  7 19  5 26  1  8  3  3 14 11 22 12 17 20 11 14 22  2 14  5 19
## [151] 22 26  8 16 25  7 23  7 16  2 22 13 19  7  4 24 11  8  9 22 25 14 18  8 22
## [176]  8 13 27 27  3  2 23  9 15 17  5 26 27 12  1  3  0 16  3  5  1  7  2  2  5

datos1[,3] # da la columna 3

##   [1] 1 2 2 1 1 1 2 2 2 1 2 2 1 1 1 1 2 2 2 2 1 2 1 1 2 1 1 1 2 2 2 1 1 1 1 2 1
##  [38] 2 2 2 2 1 1 2 2 2 1 1 2 2 1 2 1 1 1 2 2 2 1 2 2 1 1 1 1 2 2 1 1 1 2 2 2 1
##  [75] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 2 2 1 1 1 2 1 2 2 2 2 2 2 1 2 1 2 2 1 2 1
## [112] 2 1 1 2 2 1 1 2 2 1 1 1 1 1 1 2 2 2 1 1 2 2 2 2 2 1 1 1 2 1 2 2 1 1 2 1 2
## [149] 2 1 2 1 1 2 2 1 2 2 1 1 2 2 2 2 1 1 2 2 2 1 2 2 2 2 1 2 1 2 2 2 2 1 1 2 2
## [186] 2 2 1 2 2 2 2 2 1 2 2 2 1 2 2

datos1[3,] # da el renglon 3

##   antiguedad horasextra sexo cursos incapacidad aptitudes escolaridad salario
## 3         17        115    2      3           5       134           1   34875
##   edad
## 3   48

datos1$sexo  # da la columna 3

##   [1] 1 2 2 1 1 1 2 2 2 1 2 2 1 1 1 1 2 2 2 2 1 2 1 1 2 1 1 1 2 2 2 1 1 1 1 2 1
##  [38] 2 2 2 2 1 1 2 2 2 1 1 2 2 1 2 1 1 1 2 2 2 1 2 2 1 1 1 1 2 2 1 1 1 2 2 2 1
##  [75] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 2 2 1 1 1 2 1 2 2 2 2 2 2 1 2 1 2 2 1 2 1
## [112] 2 1 1 2 2 1 1 2 2 1 1 1 1 1 1 2 2 2 1 1 2 2 2 2 2 1 1 1 2 1 2 2 1 1 2 1 2
## [149] 2 1 2 1 1 2 2 1 2 2 1 1 2 2 2 2 1 1 2 2 2 1 2 2 2 2 1 2 1 2 2 2 2 1 1 2 2
## [186] 2 2 1 2 2 2 2 2 1 2 2 2 1 2 2

datos1[1:10,3:6]  # para seleccionar una parte de los datos

##    sexo cursos incapacidad aptitudes
## 1     1      4           9     121.9
## 2     2      2           2     114.2
## 3     2      3           5     134.1
## 4     1      1           1     114.0
## 5     1      2           0     151.4
## 6     1      4           3      96.7
## 7     2      2           4      98.4
## 8     2      1           1     110.1
## 9     2      2           5     102.0
## 10    1      3           2     101.0

# Si llamamos a la variable por su nombre...
#antiguedad  # habiamos llamado datos1$antiguedad

#antiguedad*12

# MARCA ERROR AL TEJER!

# Solucion.....

attach(datos1) # con el attach ya puedo accesar a los nombres de las variables

antiguedad

##   [1] 11 24 17  9 15  6  4  2 17 17 15 21  4 12 23 20 19 12  5 11 11  8 20  1  6
##  [26] 18 21  7 21 27 20 11 11  3 16  2 12 16  9 15  3 17 17 23  6  1  7  4 22 18
##  [51] 22 25 15 24 15  2 17  7  8 11 10 10  8 15  8 23 24 22  6  3 16  4 12 23 13
##  [76] 11  4  3  9 12  3 16 23  2 18  3  7 25  2 17 22  6 27 14 12 24 14 14 11  4
## [101]  3 12 12 19 12  4  5 23  9 24 21 19 14  3  3  8 15 18  7 11  5 18 12  2 26
## [126] 26 11 11  0  7 19  5 26  1  8  3  3 14 11 22 12 17 20 11 14 22  2 14  5 19
## [151] 22 26  8 16 25  7 23  7 16  2 22 13 19  7  4 24 11  8  9 22 25 14 18  8 22
## [176]  8 13 27 27  3  2 23  9 15 17  5 26 27 12  1  3  0 16  3  5  1  7  2  2  5

antiguedad*12

##   [1] 132 288 204 108 180  72  48  24 204 204 180 252  48 144 276 240 228 144
##  [19]  60 132 132  96 240  12  72 216 252  84 252 324 240 132 132  36 192  24
##  [37] 144 192 108 180  36 204 204 276  72  12  84  48 264 216 264 300 180 288
##  [55] 180  24 204  84  96 132 120 120  96 180  96 276 288 264  72  36 192  48
##  [73] 144 276 156 132  48  36 108 144  36 192 276  24 216  36  84 300  24 204
##  [91] 264  72 324 168 144 288 168 168 132  48  36 144 144 228 144  48  60 276
## [109] 108 288 252 228 168  36  36  96 180 216  84 132  60 216 144  24 312 312
## [127] 132 132   0  84 228  60 312  12  96  36  36 168 132 264 144 204 240 132
## [145] 168 264  24 168  60 228 264 312  96 192 300  84 276  84 192  24 264 156
## [163] 228  84  48 288 132  96 108 264 300 168 216  96 264  96 156 324 324  36
## [181]  24 276 108 180 204  60 312 324 144  12  36   0 192  36  60  12  84  24
## [199]  24  60

# OJO: el problema con attach es que si manejo 2 o mas bases de datos distintas,
# y estas comparten el mismo nombre de las columnas, R se hace bolas y no
# sabra exactamente a que base de datos se refiere.

#detach(datos1) # para decirle a R que ya NO usare esa base de datos

# Datos faltantes?

library(naniar)

pct_miss(datos1)

## [1] 0

1.2 Tipo de variable y escala de medición

Antes de realizar el análisis estadístico más adecuado es conveniente revisar:

1. Que la información de todas las variables sea completa
2. Que nuestra elección esté de acuerdo con el tipo de variable y su escala de medición para ¡no meter la pata!

1.2.1 Identificación de las variables

Variable	Tipo	Escala de medición
Antigüedad	cuantitativa/discreta	razón
Hora extra	cuantitativa/discreta	razon
Sexo	cualitativa	nominal
Cursos	cuantitativa/discreta	razón
Incapacidad	cuantitativa/discreta	razon
Aptitudes	cuantitativa/continua	intervalo
Escolaridad	cualitativa	ordinal
Salarios	cuantitativa/continua	razón
Edad	cuantitativa/discreta	razón

1.3 Análisis estadístico de las variables cualitativas

1.3.1 Sexo

1.3.1.1 Análisis tabular y gráfico

Basándonos en la distribución de frecuencias,

# Tabla de frecuencias

(tabla_sexo <- table(sexo))  # tabla de frecuencias absolutas

## sexo
##   1   2 
##  81 119

(tablar_sexo <- prop.table(tabla_sexo)) # tabla de frecuencias relativas

## sexo
##     1     2 
## 0.405 0.595

# otra forma de obtener las marginales:

(f.relativa_sexo<- c(sum(sexo==1)/length(sexo),sum(sexo==2)/length(sexo)))

## [1] 0.405 0.595

En resumen:

# Tabla de frecuencias completa

tabla.frec_sex<- matrix(cbind(tabla_sexo[1],tabla_sexo[2],
                              sum(tabla_sexo), tablar_sexo[1],tablar_sexo[2],sum(tablar_sexo)),
                    byrow=T,nrow = 2,ncol=3)
colnames(tabla.frec_sex)<- c("Mujer","Hombre","Totales")
rownames(tabla.frec_sex)<- c("fi","pi")
tabla.frec_sex

##     Mujer  Hombre Totales
## fi 81.000 119.000     200
## pi  0.405   0.595       1

Observamos que en esta empresa hay más empleados que empleadas. De los primeros hay 119 personas que constituyen el 59.5% del total.

Hacemos una gráfica para ver la distribución de los empleados (por género) en esta empresa:

# Diagrama de barras
barplot(tabla_sexo,names.arg=c("Mujer","Hombre"),col=c(3,32),
        ylim=c(0,140))

abline(h=tabla_sexo,col=c("magenta","blue"),lwd=2)  # dibuja las lineas horizontales
title("Diagrama de barras - Sexo \n (Frecuencias Absolutas)")

barplot(tablar_sexo,names.arg=c("M","H"),col=3:4,ylim=c(0,0.6))
title("Diagrama de barras - Sexo \n (Frecuencias Relativas)")

# Diagrama Circular

par(mfrow=c(1,2))
pie(tabla_sexo,labels=paste(c("M","H"),tabla_sexo,sep=": "),col=5:6,radius=1)
title("Diagrama de pastel - Sexo \n (Frecuencias absolutas)")   # \n para escribir en otro renglon

pie(tablar_sexo,labels=paste(c("M","H"),tablar_sexo,sep=": "),col=15:16,radius=1)
title("Diagrama de pastel - Sexo \n (Frecuencias relativas)")

De la gráfica de barras y circular se observa que, en efecto, hay diferencias entre el número de empleados y empleadas que laboran en esta empresa. Vale la pena analizar estadísticamente si esto es así.

# Comparacion entre la proporcion de mujeres y la proporcion de hombres.

n_muj<- tabla_sexo[1]
n_hom<- tabla_sexo[2]
   
prop.test(x=c(n_muj, n_hom), n=c(200,200), conf.level=0.90)

## 
##  2-sample test for equality of proportions with continuity correction
## 
## data:  c(n_muj, n_hom) out of c(200, 200)
## X-squared = 14, df = 1, p-value = 2e-04
## alternative hypothesis: two.sided
## 90 percent confidence interval:
##  -0.276 -0.104
## sample estimates:
## prop 1 prop 2 
##  0.405  0.595

De acuerdo con el intervalo de confianza para la diferencia de proporciones se concluye que SÍ hay diferencia entre el número de empleados y empleadas, de hecho \(p_{muj}-p_{hom}<0\); así que la proporción de mujeres si es estadísticamente inferior a la correspondiente de los hombres (como lo muestran las gráficas).

***NOTA: # ¿Tiene sentido la siguiente grafica?

hist(sexo,breaks=0:2)   # eje vertical es la frecuencia absoluta

hist(sexo,breaks=0:2,probability=TRUE)  # eje vertical es la frecuencia relativa

1.3.1.2 Medidas descriptivas

Sólo puede obtenerse la moda.

# install.packages(modeest)
library(modeest)    # "genefilter" es una dependencia que NO esta disponible para R 3.6.0

# Otra paqueteria
# install.packages(modes)
# library(modes)   # no disponible en R 4.0.0

# Variable nominal
mfv(sexo)     # most frequent value

## [1] 2

mode(sexo)    # NO!

## [1] "numeric"

# modes(sex)  # library(modes) 

summary(sexo)   # no tiene interpretacion

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1.00    1.00    2.00    1.59    2.00    2.00

1.3.2 Escolaridad

1.3.2.1 Análisis tabular y gráfico

Basándonos en la distribución de frecuencias,

# Tabla de frecuencias

(tabla_esc <- table(escolaridad))  # tabla de frecuencias absolutas

## escolaridad
##   0   1   2   3 
##  34 117  47   2

(tablar_esc <- prop.table(tabla_esc)) # tabla de frecuencias relativas

## escolaridad
##     0     1     2     3 
## 0.170 0.585 0.235 0.010

En resumen:

# Tabla de frecuencias completa

tabla.frec_esc<- matrix(cbind(tabla_esc[1],tabla_esc[2],
                              tabla_esc[3], tabla_esc[4],
                              sum(tabla_esc),
                              tablar_esc[1],tablar_esc[2],
                              tablar_esc[3], 
                              tablar_esc[4], sum(tablar_esc)),
                    byrow=T,nrow = 2,ncol=5)
colnames(tabla.frec_esc)<- c("B","L(st)","L(ct)","P","Totales")
rownames(tabla.frec_esc)<- c("fi","pi")
tabla.frec_esc

##        B   L(st)  L(ct)    P Totales
## fi 34.00 117.000 47.000 2.00     200
## pi  0.17   0.585  0.235 0.01       1

Observamos que en esta empresa hay 117 emplead@s con licenciatura SIN título y representan el 58.5% del total. La diferencia entre los que sólo tienen el bachillerato y licenciatura CON título es pequeña (y veremos sí es significativa). Definitivamente casi no hay emplead@s con posgrado.

Hacemos una gráfica para ver la distribución de los empleados (por escolaridad) en esta empresa:

# Diagrama de barras
barplot(tabla_esc,names.arg=c("B","L(st)","L(ct)","P"),
        col=c(3,32,2),ylim=c(0,140))

# Diagrama Circular

pie(tablar_esc,labels=paste(c("B","L(st)","L(ct)","P"),
                            tablar_esc,
                            sep=": "),col=rainbow(4),radius=1)
title("Diagrama de pastel - Escolaridad \n (Frecuencias relativas)")

De la gráfica de barras y circular se observa que, en efecto, hay diferencias entre el nivel educativo de los empleado(a)s. Vale la pena analizar estadísticamente si esto es así.

# Comparacion entre la proporcion de mujeres y la proporcion de hombres.

n_bach<- tabla_esc[1]
n_lct<- tabla_esc[3]
   
prop.test(x=c(n_bach, n_lct), n=c(200,200), conf.level=0.90)

## 
##  2-sample test for equality of proportions with continuity correction
## 
## data:  c(n_bach, n_lct) out of c(200, 200)
## X-squared = 2, df = 1, p-value = 0.1
## alternative hypothesis: two.sided
## 90 percent confidence interval:
##  -0.13588  0.00588
## sample estimates:
## prop 1 prop 2 
##  0.170  0.235

Basándonos en el valor-p, NO se rechaza H0 por lo que \(p_{bach}-p_{lct}=0\) y por consiguiente no hay evidencia estadística suficiente para decir que la proporción de emplead@s con licenciatura CON título es diferente a la correspondiente con bachillerato.

1.3.2.2 Medidas descriptivas

Sólo tienen sentido: moda y mediana (o cualquier porcentil)

1. Moda:

mfv(escolaridad)

## [1] 1

Los emplead@s con licenciatura SIN titulo son los que más laboran en esta empresa

2. Mediana

median(escolaridad)

## [1] 1

quantile(escolaridad,prob=c(0.10,0.25,0.5,0.75))

## 10% 25% 50% 75% 
##   0   1   1   1

El 10% de los emplead@s tienen a lo más bachillerato. Así que el 90% restante tienen educación superior. WOW!

El 75% de los emplead@s tienen a lo más licenciatura SIN título, por lo que el 25% restante ya tienen el título o poseen un posgrado.

1.4 Análisis estadístico de las variables cuantitativas

1.4.1 Cursos

1.4.1.1 Análisis tabular y gráfico

Se trata de una variable discreta. Así que podemos presentar su resumen mediante una distribución de frecuencias como la correspondiente a una variable cualitativa.

# Tabla de frecuencias

(tabla_cursos <- table(cursos))  # tabla de frecuencias absolutas

## cursos
##  0  1  2  3  4  5  6  9 
## 27 41 43 41 27 19  1  1

(tablar_cursos <- prop.table(tabla_cursos)) # tabla de frecuencias relativas

## cursos
##     0     1     2     3     4     5     6     9 
## 0.135 0.205 0.215 0.205 0.135 0.095 0.005 0.005

En resumen

# Tabla de frecuencias completa

tabla.frec_cursos<- matrix(cbind(tabla_cursos[1],tabla_cursos[2],
                                tabla_cursos[3], tabla_cursos[4],
                                tabla_cursos[5], tabla_cursos[6],
                                tabla_cursos[7], tabla_cursos[8],
                                sum(tabla_cursos),
                                tablar_cursos[1], tablar_cursos[2],
                                tablar_cursos[3], tablar_cursos[4],
                                tablar_cursos[5], tablar_cursos[6],
                                tablar_cursos[7], tablar_cursos[8],
                                sum(tablar_cursos)),
                                byrow=T,nrow = 2,ncol=9)
colnames(tabla.frec_cursos)<- c("0","1","2","3","4","5","6","9","Totales")
rownames(tabla.frec_cursos)<- c("fi","pi")
tabla.frec_cursos

##         0      1      2      3      4      5     6     9 Totales
## fi 27.000 41.000 43.000 41.000 27.000 19.000 1.000 1.000     200
## pi  0.135  0.205  0.215  0.205  0.135  0.095 0.005 0.005       1

Veamos la distribución del número de cursos para los emplead@s de esta empresa.

par(mfrow=c(1,2))
barplot(tabla_cursos)  # frecuencias absolutas
title("Diagrama de barras (Cursos) \n (Frecuencias absolutas)")

barplot(tablar_cursos,col=c(10:20))   # frecuencias relativas
title("Diagrama de barras (Points) \n (Frecuencias relativas)")

¿Qúe pasa si a la ingenua graficamos un histograma?

hist(cursos,breaks=-1:9,plot=FALSE)  # las clases son: [-1,0] (0,1],...,(8,9])

## $breaks
##  [1] -1  0  1  2  3  4  5  6  7  8  9
## 
## $counts
##  [1] 27 41 43 41 27 19  1  0  0  1
## 
## $density
##  [1] 0.135 0.205 0.215 0.205 0.135 0.095 0.005 0.000 0.000 0.005
## 
## $mids
##  [1] -0.5  0.5  1.5  2.5  3.5  4.5  5.5  6.5  7.5  8.5
## 
## $xname
## [1] "cursos"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"

# OJO: las clases (breaks) los defini asi para que se pareciera al diagrama
# de barras y poderlo comparar.

hist(cursos,breaks=-1:9,probability=TRUE)

El perfil en ambas gráficas es “idéntico”, sólo que en el segundo caso las barras están pegadas. ¿Será correcto?

Podemos pensar que los valores de esta variable discreta son como las categorías y, por consiguiente, el diagrama de pay sería adecuado:

par(mfrow=c(1,2))
pie(tabla_cursos,labels=paste(c("0","1","2","3","4","5","6","9"),
                              tabla_cursos,sep=": "),col=palette("Alphabet"),radius=1)
title("Diagrama de pay - Cursos \n (Frecuencias absolutas)")

pie(tablar_cursos,labels=paste(c("0","1","2","3","4","5","6","9"),
                               tablar_cursos,sep=": "),col=topo.colors(8),radius=1)
title("Diagrama de pay - Cursos \n (Frecuencias relativas)")

Por tratarse de una variable cuantitativa es posible representar sus valores mediante un diagrama de puntos (dotplot o strip chart).

dotchart(cursos,labels=row.names(datos1),cex=0.8,groups = as.factor(sexo),#forzar a sexo a ser categorias, si son, pero no lo sabe R
         xlab="Número de cursos",main="Diagrama de puntos")  # Cleveland dot plot

Mejor tratemos con un subconjunto:

datos2<- subset(datos1,cursos>4); datos2

##     antiguedad horasextra sexo cursos incapacidad aptitudes escolaridad salario
## 32          11        140    1      5           7      94.8           0   27500
## 41           3          6    2      9           3     107.9           2   19000
## 46           1         20    2      5           4     104.7           2   27050
## 49          22         77    2      5           4      80.5           0   40900
## 66          23         74    2      5           4     110.2           1   40370
## 67          24        138    2      5           4     103.6           1   42450
## 73          12          0    2      6           8     104.8           1   35850
## 76          11         90    2      5           2     101.7           2   40700
## 80          12        136    2      5           4     113.2           2   41800
## 87           7        132    2      5           7     126.3           1   28215
## 118         18         44    1      5           5      96.2           1   36200
## 122         18        182    1      5           8     104.4           1   36100
## 131         19         13    1      5           3     119.1           2   38300
## 143         20         39    2      5           2     123.4           1   43945
## 147          2         90    1      5           6     103.5           1   23325
## 153          8        169    1      5           4      87.8           2   31625
## 166         24         37    1      5           4      99.1           1   47125
## 169          9        198    2      5           3     107.8           1   31570
## 177         13        129    1      5           4     104.2           1   38300
## 180          3        209    2      5           3      92.4           0   23600
## 190          1         70    2      5           4      96.4           0   20000
##     edad
## 32    32
## 41    23
## 46    31
## 49    63
## 66    65
## 67    63
## 73    61
## 76    42
## 80    42
## 87    31
## 118   51
## 122   46
## 131   40
## 143   62
## 147   25
## 153   34
## 166   64
## 169   44
## 177   38
## 180   22
## 190   20

dotchart(datos2$cursos,labels=row.names(datos2),cex=0.8,xlab="Número de cursos",main="Diagrama de puntos")  # Cleveland dot plot

¿Y si mejor seleccionamos al azar?

#si queremos que sea lo mismo poner una semilla
datos3<- sample(datos1$cursos,size=60,replace = F); datos3 #replace F  es igual a sin remplazo

##  [1] 3 9 3 1 1 2 1 1 1 1 3 3 4 1 4 0 0 5 2 4 1 1 2 2 3 3 3 1 2 0 3 0 1 4 6 1 4 3
## [39] 1 3 5 2 4 2 5 2 2 1 2 1 1 1 5 2 2 3 3 4 4 5

dotchart(datos3,cex=0.8,xlab="Numero de cursos",main="Diagrama de puntos")  # Cleveland dot plot

# Muestreo aleatorio de TODO el conjunto de datos

#sample(datos1,10,replace = F)  # no es correcto!

datos4<- datos1[sample(nrow(datos1),30,replace = F),]; datos4  

##     antiguedad horasextra sexo cursos incapacidad aptitudes escolaridad salario
## 164          7        123    2      4           2     111.0           1   31575
## 56           2        169    2      1           9     118.4           1   15625
## 20          11         27    2      2           4     101.7           1   23960
## 31          20         99    2      0           4      96.9           2   37100
## 129          0         81    2      3           7     104.2           1   33900
## 30          27        157    2      2           5     113.2           2   39975
## 148         14         37    2      4           3     106.4           2   31325
## 106          4        200    1      2           3     100.9           2   32350
## 39           9        124    2      1           5      87.7           2   27200
## 168          8        100    2      2           6     107.8           2   40000
## 179         27        118    2      0           8     119.1           2   44000
## 149          5        173    2      2           4     110.2           1   34325
## 84           2        187    2      3           4     121.5           2   25150
## 122         18        182    1      5           8     104.4           1   36100
## 200          5         74    2      3           6     111.0           1   33000
## 69           6         68    1      2           5      87.6           1   39500
## 170         22        198    1      1           9     122.7           0   41750
## 189         12        209    2      2           3     117.5           1   31450
## 196          1         11    2      1           4      93.9           3   30200
## 159         16        105    1      2           5      89.7           1   39400
## 6            6         43    1      4           3      96.7           1   22635
## 90          17        118    2      1           3      98.1           0   40425
## 54          24        215    1      0           7      95.3           0   26180
## 124          2        217    1      2           0     123.1           1   34550
## 153          8        169    1      5           4      87.8           2   31625
## 107          5        182    2      1           7     124.1           1   34900
## 142         17         28    2      0           3     105.1           2   37675
## 94          14        193    1      3           4     122.0           1   42450
## 150         19          6    1      2           4     105.0           1   40500
## 123         12         48    1      4           3     111.4           1   28150
##     edad
## 164   34
## 56    26
## 20    30
## 31    57
## 129   31
## 30    67
## 148   32
## 106   31
## 39    27
## 168   41
## 179   66
## 149   35
## 84    29
## 122   46
## 200   34
## 69    57
## 170   40
## 189   34
## 196   30
## 159   41
## 6     45
## 90    45
## 54    41
## 124   31
## 153   34
## 107   35
## 142   42
## 94    57
## 150   50
## 123   35

head(datos4)

##     antiguedad horasextra sexo cursos incapacidad aptitudes escolaridad salario
## 164          7        123    2      4           2     111.0           1   31575
## 56           2        169    2      1           9     118.4           1   15625
## 20          11         27    2      2           4     101.7           1   23960
## 31          20         99    2      0           4      96.9           2   37100
## 129          0         81    2      3           7     104.2           1   33900
## 30          27        157    2      2           5     113.2           2   39975
##     edad
## 164   34
## 56    26
## 20    30
## 31    57
## 129   31
## 30    67

dotchart(datos4$cursos,labels=row.names(datos4),cex=0.8,xlab="Numero de cursos",main="Diagrama de puntos")  # Cleveland dot plot

Otro tipo de gráfica de puntos conocida como diagrama de puntos

#install.packages("BHH2")
library(BHH2)

dotPlot(na.omit(datos4$cursos),main="Diagrama de puntos (muestra aleatoria)",xlab="Número de puntos",pch=16)

Analicemos el diagrama de tallo y hojas. Recordemos que en este diagrama se divide a cada observación en dos partes: tallo y hoja.

1. Cuando hay punto decimal, éste se usa como separador.
2. El problema es cuando se tienen muchos dígitos a la derecha del separador (hojas). Esto es poco práctico, por lo que sólo se considera un decimal, ¿cómo escogemos este único decimal?

El argumento “scale” permite modificar la longitud del diagrama y con ello tener este único decimal. Por ejemplo:

stem(datos1$cursos, scale=0.1)  

## 
##   The decimal point is 1 digit(s) to the right of the |
## 
##   0 | 00000000000000000000000000011111111111111111111111111111111111111111+99
##   0 | 555555555555555555569

Los valores que se leen son: 0.0,0.1,0.5,0.6,0.9

stem(datos1$cursos, scale=1) #no se usa mucho en ciendia de datos

## 
##   The decimal point is at the |
## 
##   0 | 000000000000000000000000000
##   0 | 
##   1 | 00000000000000000000000000000000000000000
##   1 | 
##   2 | 0000000000000000000000000000000000000000000
##   2 | 
##   3 | 00000000000000000000000000000000000000000
##   3 | 
##   4 | 000000000000000000000000000
##   4 | 
##   5 | 0000000000000000000
##   5 | 
##   6 | 0
##   6 | 
##   7 | 
##   7 | 
##   8 | 
##   8 | 
##   9 | 0

Los valores que se leen son: 0,1,2,3,4,5,6,9

COMENTARIO: En el primer caso, el valor exacto se pierde. Sin embargo, esto no es relevante para este tipo de gráfica/tabla. Sólo interesa identificar alrededor de cual valor se concentran más las observaciones y el perfil que muestra esta distribución.

Finalmente analizamos el diagrama de caja y brazos:

boxplot(cursos,horizontal=T,col="yellow") 

Se observa que hay una persona con más de 8 cursos y es un outlier. Claramente la distribución de cursos es asimétrica hacia la derecha. El 50% de los emplead@s ha tomado a lo más 2 cursos; mientras que el 75% ha llevado hasta 3 cursos. Veamos de quién se trata para premiarlo.

datos1[datos1$cursos>8,]

##    antiguedad horasextra sexo cursos incapacidad aptitudes escolaridad salario
## 41          3          6    2      9           3       108           2   19000
##    edad
## 41   23

datos1[datos1$cursos>8,]$salario

## [1] 19000

1.4.1.2 Medidas descriptivas

1. Moda

library(fdth)

## 
## Attaching package: 'fdth'

## The following object is masked from 'package:modeest':
## 
##     mfv

## The following objects are masked from 'package:stats':
## 
##     sd, var

mfv(cursos) #Usa libreria obligatoriamente

## [1] 2

2. Media

mean(cursos)

## [1] 2.34

3. Cuantiles y Amplitud intercuartílica

quantile(cursos, probs = c(0.25,0.50,0.75))

## 25% 50% 75% 
##   1   2   3

IQR(cursos)

## [1] 2

4. Varianza

var(cursos); sd(cursos)

## [1] 2.57

## [1] 1.6

5. Error medio

mean(abs(na.omit(cursos)-mean(na.omit(cursos))))   # error medio (media)

## [1] 1.33

mean(abs(cursos-median(cursos)))  # error medio (mediana)

## [1] 1.29

mean(abs(cursos-mfv(cursos)))  # error medio (moda) OJO!!!

## [1] 1.29

6. Coeficiente de variación

c.v<- sd(cursos)/mean(cursos); c.v # tenemos poca variavilidad

## [1] 0.685

7. Asimetría y curtosis

A veces es difícil identificar si la distribución de frecuencias es asimétrica o simétrica mediante un histograma o un diagrama de barras. Para ello existen medidas estadísticas que nos permiten saber si el perfil es simétrico o no; además si es un perfil “picudo” o “plano”. Estas medidas son: coeficiente de asimetría y coeficiente de curtosis.

1. Coeficiente de asimetría

n<- length(cursos)   # longitud del vector
mu3<- sum((cursos-mean(cursos))^3)/(n-1)
s3<- sd(cursos)^3
j3<- mu3/s3; j3    

## [1] 0.478

# Usando R

library(agricolae) #se puede con libreria momentos

## 
## Attaching package: 'agricolae'

## The following object is masked from 'package:modeest':
## 
##     skewness

(j3_a<- skewness(cursos)) 

## [1] 0.483

(j3_b<- moments::skewness(cursos))

## [1] 0.479

2. Coeficiente de curtosis

mu4<- sum((cursos-mean(cursos))^4)/(n-1)
s4<- sd(cursos)^4
j4<- mu4/s4; j4    

## [1] 3.24

# Usando R

library(agricolae) #numero 3 distribucion normal
(j4_a<- kurtosis(cursos)+3)  # da el exceso de curtosis  

## [1] 3.3

(j4_b<- moments::kurtosis(aptitudes))  

## [1] 3.39

O bien, podemos resumir parte de esta información en:

fivenum(cursos) # minimum, lower-hinge, median, upper-hinge, maximum)

## [1] 0 1 2 3 9

summary(cursos)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    1.00    2.00    2.34    3.00    9.00

---

---

COMENTARIO SOBRE “HINGES”: es similar a los cuartiles. En conjuntos grandes de datos, la diferencia entre ellos es muy pequeña. Según la definición de Tukey los “hinges” son:

1. Hinges inferior: es la mediana de la mitad inferior de todos los valores a la izquierda de la mediana de los datos.

2. Hinges superior: es la mediana de la mitad superior de todos los valores a la derecha de la mediana de los datos

Ejemplo 1: n=20 observaciones La mediana está en la posición 10.5. La mitad inferior comprende las posiciones 1 hasta 10; mientras que la superior desde la 11 hasta 20.

El “lower hinge” es la mediana de la mitad inferior y sería la posición 5.5.

El “upper hinge” es la mediana de la mitad superior y sería la posición 5.5 (comenzando desde la posición 11), por lo que a final de cuentas será la posición 15.5.

Ejemplo 2: n=21 observaciones

La mediana está en la posición 11. La mitad inferior comprende las posiciones 1 hasta 11; mientras que la superior desde la 11 hasta 21.

El “lower hinge” es la mediana de la mitad inferior y sería la posición 6.

El “upper hinge” es la mediana de la mitad superior y sería la posición 6 (comenzando desde la posición 11), por lo que a final de cuentas será la posición 16.

---

---

1.4.2 Salario

1.4.2.1 Análisis tabular y gráfico

Se trata de una variable continua, así que si hacemos lo siguiente:

table(salario)

## salario
## 15625 18955 19000 19575 19775 20000 20430 22000 22025 22400 22635 23065 23100 
##     1     1     1     1     1     2     1     1     1     1     1     1     2 
## 23150 23325 23550 23600 23650 23685 23960 24450 24525 24650 24750 24850 25100 
##     1     1     1     3     1     1     1     1     1     2     1     1     2 
## 25150 25545 25595 26180 26250 26275 26500 26800 27050 27175 27180 27200 27300 
##     1     1     1     1     1     1     3     1     1     1     1     1     1 
## 27350 27500 27600 27940 27950 28090 28150 28215 28650 28750 29000 29100 29500 
##     1     3     1     1     1     1     2     1     1     1     2     1     1 
## 29900 30000 30200 30440 30495 30500 30600 30630 31200 31325 31400 31450 31570 
##     1     1     1     1     1     1     1     1     1     1     1     2     1 
## 31575 31625 32000 32350 32670 32775 33000 33110 33550 33900 33965 33975 34250 
##     1     1     3     1     1     1     3     1     1     1     1     1     1 
## 34325 34450 34550 34850 34875 34900 34975 35000 35200 35300 35600 35850 35875 
##     2     1     2     1     2     1     1     2     2     1     1     2     1 
## 36100 36195 36200 36475 37000 37100 37200 37225 37300 37375 37600 37675 37750 
##     2     1     1     1     1     1     1     1     1     1     2     1     2 
## 38170 38300 38500 38790 38900 38950 39275 39375 39400 39500 39750 39975 40000 
##     1     3     2     1     1     1     1     1     1     3     1     1     1 
## 40050 40370 40400 40425 40500 40700 40900 41150 41575 41600 41650 41705 41750 
##     1     1     1     1     3     1     1     1     1     2     1     1     1 
## 41800 41900 41910 42000 42250 42450 42850 43175 43650 43675 43780 43800 43850 
##     3     1     1     1     1     2     1     1     2     1     1     2     2 
## 43945 44000 45000 45100 45375 45800 46100 46550 46900 47025 47100 47125 47200 
##     1     2     1     1     2     1     1     1     1     1     1     1     1 
## 47300 47950 
##     1     1

barplot(salario)

COMENTARIO: Se observa que no es util para este tipo de variable representarlo asi, por eso conviene hacer una tabla de frecuencias usando INTERVALOS DE CLASE.

En tenemos dos opciones para construir una tabla de frecuencias para datos agrupados.

1. Usando la funcion “hist”.

Sugerencia: revisar primero los valores mínimo y máximo

#imagen del logo en sea, guardado en el directorio de trabajo
range(salario)   # minimo y maximo

## [1] 15625 47950

Con ellos, ahora sí definimos los breaks del histograma:

hist_sal <- hist(datos1$salario,breaks=seq(15600,49000,2500),plot=F)

summary(hist_sal)   

##          Length Class  Mode     
## breaks   14     -none- numeric  
## counts   13     -none- numeric  
## density  13     -none- numeric  
## mids     13     -none- numeric  
## xname     1     -none- character
## equidist  1     -none- logical

Construimos la correspondiente tabla de frecuencias

n <- length(hist_sal$breaks)
tab_sal <- cbind(hist_sal$breaks[-n],hist_sal$breaks[-1],
                 hist_sal$counts,
                 hist_sal$counts/sum(hist_sal$counts),
                 cumsum(hist_sal$counts),
                 cumsum(hist_sal$counts/sum(hist_sal$counts)))
dimnames(tab_sal)[[2]]<-c("Linf","Lsup","f","fr","F","Fr")
print(tab_sal)

##        Linf  Lsup  f    fr   F    Fr
##  [1,] 15600 18100  1 0.005   1 0.005
##  [2,] 18100 20600  7 0.035   8 0.040
##  [3,] 20600 23100  7 0.035  15 0.075
##  [4,] 23100 25600 20 0.100  35 0.175
##  [5,] 25600 28100 20 0.100  55 0.275
##  [6,] 28100 30600 16 0.080  71 0.355
##  [7,] 30600 33100 18 0.090  89 0.445
##  [8,] 33100 35600 22 0.110 111 0.555
##  [9,] 35600 38100 19 0.095 130 0.650
## [10,] 38100 40600 25 0.125 155 0.775
## [11,] 40600 43100 19 0.095 174 0.870
## [12,] 43100 45600 16 0.080 190 0.950
## [13,] 45600 48100 10 0.050 200 1.000

2. Usando el paquete “fdth”

library(fdth)
(tabla_cont<- fdt(salario,breaks="Scott",start=15600,end=49000,h=2500,right=T))

##   Class limits  f   rf rf(%)  cf cf(%)
##  (15600,18100]  1 0.00   0.5   1   0.5
##  (18100,20600]  7 0.04   3.5   8   4.0
##  (20600,23100]  7 0.04   3.5  15   7.5
##  (23100,25600] 20 0.10  10.0  35  17.5
##  (25600,28100] 20 0.10  10.0  55  27.5
##  (28100,30600] 16 0.08   8.0  71  35.5
##  (30600,33100] 18 0.09   9.0  89  44.5
##  (33100,35600] 22 0.11  11.0 111  55.5
##  (35600,38100] 19 0.10   9.5 130  65.0
##  (38100,40600] 25 0.12  12.5 155  77.5
##  (40600,43100] 19 0.10   9.5 174  87.0
##  (43100,45600] 16 0.08   8.0 190  95.0
##  (45600,48100] 10 0.05   5.0 200 100.0

(tabla_cont<- fdt(salario,breaks="Sturges",right=T))

##         Class limits  f   rf rf(%)  cf cf(%)
##  (15468.75,19131.06]  3 0.01   1.5   3   1.5
##  (19131.06,22793.36]  9 0.04   4.5  12   6.0
##  (22793.36,26455.67] 26 0.13  13.0  38  19.0
##  (26455.67,30117.97] 28 0.14  14.0  66  33.0
##  (30117.97,33780.28] 25 0.12  12.5  91  45.5
##  (33780.28,37442.58] 34 0.17  17.0 125  62.5
##  (37442.58,41104.89] 32 0.16  16.0 157  78.5
##  (41104.89,44767.19] 29 0.14  14.5 186  93.0
##   (44767.19,48429.5] 14 0.07   7.0 200 100.0

(tabla_cont<- fdt(salario,breaks="Sturges",start=15600,end=49000,h=2500,right=T))

##   Class limits  f   rf rf(%)  cf cf(%)
##  (15600,18100]  1 0.00   0.5   1   0.5
##  (18100,20600]  7 0.04   3.5   8   4.0
##  (20600,23100]  7 0.04   3.5  15   7.5
##  (23100,25600] 20 0.10  10.0  35  17.5
##  (25600,28100] 20 0.10  10.0  55  27.5
##  (28100,30600] 16 0.08   8.0  71  35.5
##  (30600,33100] 18 0.09   9.0  89  44.5
##  (33100,35600] 22 0.11  11.0 111  55.5
##  (35600,38100] 19 0.10   9.5 130  65.0
##  (38100,40600] 25 0.12  12.5 155  77.5
##  (40600,43100] 19 0.10   9.5 174  87.0
##  (43100,45600] 16 0.08   8.0 190  95.0
##  (45600,48100] 10 0.05   5.0 200 100.0

(tabla_cont<- fdt(salario,breaks="FD",right=T))

##         Class limits  f   rf rf(%)  cf cf(%)
##  (15468.75,19588.84]  4 0.02   2.0   4   2.0
##  (19588.84,23708.94] 19 0.10   9.5  23  11.5
##  (23708.94,27829.03] 29 0.14  14.5  52  26.0
##  (27829.03,31949.12] 28 0.14  14.0  80  40.0
##  (31949.12,36069.22] 34 0.17  17.0 114  57.0
##  (36069.22,40189.31] 35 0.17  17.5 149  74.5
##  (40189.31,44309.41] 37 0.18  18.5 186  93.0
##   (44309.41,48429.5] 14 0.07   7.0 200 100.0

Hacemos el histograma:

par(mfrow=c(1,2))
hist(salario,col=rainbow(6), prob=T)
hist(salario,col=rainbow(6), prob=F)

El histograma muestra el perfil de la distribución de probabilidad que mejor describe a la variable. Para ello, resulta conveniente dibujar el polígono de frecuencias. De hecho,

hist(salario,probability=T,col="grey50")  # histograma cuyo eje vertical esta en escala
# de "probabilidades".

lines(density(salario),lwd=2,col="blue")   # curva de densidad de probabilidad

par(mfrow=c(1,2))
library(MASS)
frequency.polygon(salario,nclass = nclass.freq(salario))

## Warning in hist.default(x, nclass, probability = TRUE, plot = FALSE, ...):
## argument 'probability' is not made use of

library(UsingR)

## Loading required package: HistData

## Loading required package: Hmisc

## Loading required package: lattice

## Loading required package: survival

## Loading required package: Formula

## Loading required package: ggplot2

## 
## Attaching package: 'Hmisc'

## The following objects are masked from 'package:base':
## 
##     format.pval, units

## 
## Attaching package: 'UsingR'

## The following object is masked from 'package:survival':
## 
##     cancer

simple.freqpoly(salario,col=heat.colors(6))

NOTA: si se requieren intervalos desiguales para observar ciertos detalles de interés

hist_uneq<- hist(salario,breaks = quantile(salario, 0:10 / 10),col="gray")  # usa los cuantiles 

hist(salario,breaks=c(15600,seq(27000,32000,800),37690,49000))   # intervalos en puntos arbitrarios

Aprovechando la tabla de frecuencias con intervalos de clase y su histograma, vale la pena revisar la siguiente gráfica: ojiva

library(agricolae)
ogive.freq(hist_sal,main="Ojiva (salario)",xlab="Límite superior de
           la clase",ylab="Pi (%)")   

##    Límite superior de\n           la clase   RCF
## 1                                    15600 0.000
## 2                                    18100 0.005
## 3                                    20600 0.040
## 4                                    23100 0.075
## 5                                    25600 0.175
## 6                                    28100 0.275
## 7                                    30600 0.355
## 8                                    33100 0.445
## 9                                    35600 0.555
## 10                                   38100 0.650
## 11                                   40600 0.775
## 12                                   43100 0.870
## 13                                   45600 0.950
## 14                                   48100 1.000
## 15                                   50600 1.000

# "ogive.freq" requiere al histograma como objeto 

¿Cómo se interpreta esta gráfica?

# Haremos una grafica interactiva para Markdown. 
# Primero debemos crear un data.frame con las columnas de la tabla de # frecuencias (requisito de la funcion "plot_ly") 

library(plotly)

## Registered S3 method overwritten by 'httr':
##   method         from  
##   print.response rmutil

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:Hmisc':
## 
##     subplot

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:MASS':
## 
##     select

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

h<- data.frame(cbind(tab_sal[,2],tab_sal[,6])); h

##       X1    X2
## 1  18100 0.005
## 2  20600 0.040
## 3  23100 0.075
## 4  25600 0.175
## 5  28100 0.275
## 6  30600 0.355
## 7  33100 0.445
## 8  35600 0.555
## 9  38100 0.650
## 10 40600 0.775
## 11 43100 0.870
## 12 45600 0.950
## 13 48100 1.000

plot_ly(data=h,x= ~tab_sal[,2],y= ~tab_sal[,6], 
        marker = list(size = 10, color="red"),
        type = 'scatter', mode = 'lines')%>%
  layout(title = 'Ojiva (Salario)', 
         xaxis = list(title = 'Límite superior'), 
         yaxis = list(title = 'Pi'),
         shapes = list(
                  # linea vertical 
                  list(type = "line", x0 = 0, x1 = 0, 
                      y0 = 0, y1 = 1, yref = "paper"),
                  # linea horizontal 
                    list(type = "line", x0 = 0, x1 = 50000, 
                            y0 = 1, y1 = 1, xref = "paper")))

## A marker object has been specified, but markers is not in the mode
## Adding markers to the mode...

Veamos otras gráficas que ya habíamos realizado en el caso discreto.

dotPlot(salario,pch=16,col=5)   # libreria(BHH2)

stem(salario)

## 
##   The decimal point is 3 digit(s) to the right of the |
## 
##   14 | 6
##   16 | 
##   18 | 0068
##   20 | 004
##   22 | 004611123666677
##   24 | 055778911256
##   26 | 233555812223455569
##   28 | 012227800159
##   30 | 024556623455666
##   32 | 000478000169
##   34 | 0033356699990002236999
##   36 | 1122501223466788
##   38 | 2333558903445558
##   40 | 00144455579266677888899
##   42 | 035592777888999
##   44 | 0001448
##   46 | 16901123
##   48 | 0

boxplot(salario,horizontal = T,col="lightgreen")

# Mejor...

library(plotly)
plot_ly(datos1,x = ~salario, type = "box")

1.4.2.2 Medidas descriptivas

1. Moda

¿Qué resultado muestra lo siguiente?

mfv(salario)

## [1] 27500

¿Tiene sentido para una variable continua?

Nótese que si las observaciones se redondean es posible que encontremos un valor que se repita varias veces. El programa R lo que hace es redondear hacia abajo (floor) y con ello localiza la moda. Esto no es adecuado.

Si conocemos la distribución de probabilidad para la variable “salario”, entonces basta con derivarla para tener el valor que maximiza su probabilidad. ¿Y si no, qué hacemos?

Respuesta: usamos su tabla de frecuencias

# Habiamos visto que

tabla_cont   # tabla de frecuencias para datos agrupados por intervalos de clase

##         Class limits  f   rf rf(%)  cf cf(%)
##  (15468.75,19588.84]  4 0.02   2.0   4   2.0
##  (19588.84,23708.94] 19 0.10   9.5  23  11.5
##  (23708.94,27829.03] 29 0.14  14.5  52  26.0
##  (27829.03,31949.12] 28 0.14  14.0  80  40.0
##  (31949.12,36069.22] 34 0.17  17.0 114  57.0
##  (36069.22,40189.31] 35 0.17  17.5 149  74.5
##  (40189.31,44309.41] 37 0.18  18.5 186  93.0
##   (44309.41,48429.5] 14 0.07   7.0 200 100.0

str(tabla_cont)    # para ver como debemos llamar a cada elemento

## List of 2
##  $ table :'data.frame':  8 obs. of  6 variables:
##   ..$ Class limits: Factor w/ 8 levels "(15468.75,19588.84]",..: 1 2 3 4 5 6 7 8
##   ..$ f           : int [1:8] 4 19 29 28 34 35 37 14
##   ..$ rf          : num [1:8] 0.02 0.095 0.145 0.14 0.17 0.175 0.185 0.07
##   ..$ rf(%)       : num [1:8] 2 9.5 14.5 14 17 17.5 18.5 7
##   ..$ cf          : num [1:8] 4 23 52 80 114 149 186 200
##   ..$ cf(%)       : num [1:8] 2 11.5 26 40 57 74.5 93 100
##  $ breaks: Named num [1:4] 15469 48430 4120 1
##   ..- attr(*, "names")= chr [1:4] "start" "end" "h" "right"
##  - attr(*, "class")= chr [1:3] "fdt.default" "fdt" "list"

which.max(tabla_cont$table$f)  # indica el renglon cuya frecuencia absoluta, f, es maxima

## [1] 7

tabla_cont$table$`Class limits`[which.max(tabla_cont$table$f)] # indica la CLASE MODAL

## [1] (40189.31,44309.41]
## 8 Levels: (15468.75,19588.84] (19588.84,23708.94] ... (44309.41,48429.5]

Comparamos los resultados y claramente son ¡MUY DISTINTOS!

2. Media

mean(salario)

## [1] 34009

3. Cuantiles y Amplitud intercuartílica

quantile(salario, probs = c(0.25,0.50,0.75))

##   25%   50%   75% 
## 27500 34700 40378

IQR(salario)

## [1] 12878

4. Varianza

var(salario); sd(salario)

## [1] 57075260

## [1] 7555

5. Error medio

mean(abs(na.omit(salario)-mean(na.omit(salario))))   # error medio (media)

## [1] 6422

mean(abs(salario-median(salario)))  # error medio (mediana)

## [1] 6398

mean(abs(salario-mfv(salario)))  # error medio (moda) OJO!!!

## [1] 8245

6. Coeficiente de variación

c.v<- sd(salario)/mean(salario); c.v

## [1] 0.222

7. Asimetría y curtosis

1. Coeficiente de asimetría

n<- length(salario)   # longitud del vector
mu3<- sum((salario-mean(salario))^3)/(n-1)
s3<- sd(salario)^3
j3<- mu3/s3; j3    

## [1] -0.143

# Usando R

library(agricolae)
(j3_a<- skewness(salario)) 

## [1] -0.144

(j3_b<- moments::skewness(salario))

## [1] -0.143

2. Coeficiente de curtosis

mu4<- sum((salario-mean(salario))^4)/(n-1)
s4<- sd(salario)^4
j4<- mu4/s4; j4    

## [1] 2.03

# Usando R

library(agricolae)
(j4_a<- kurtosis(salario)+3)  # da el exceso de curtosis  

## [1] 2.05

(j4_b<- moments::kurtosis(salario))  

## [1] 2.04

O bien, podemos resumir parte de esta información en:

fivenum(salario) # minimum, lower-hinge, median, upper-hinge, maximum)

## [1] 15625 27500 34700 40385 47950

summary(salario)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   15625   27500   34700   34009   40378   47950

1.5 Comparación y asociación entre variables

1.5.1 Variable nominal/ordinal vs variable cuantitativa: COMPARACION ENTRE VARIABLES

Para comparar el comportamiento de estas variables conviene usar un diagrama de caja-brazos. Este diagrama se usa para la variable cuantitativa, la cualitativa solo sirve para separar el análisis por casos.

Escogemos las variables escolaridad y salario.

El diagrama de caja-brazos de salario se muestra en la sección 2.4.2.1. Es claro que hay asimetría hacia la izquierda (\(\gamma_{3}=-0.143\)) pero con esta gráfica aún no podemos saber si la escolaridad juega un papel relevante. Para ello hacemos el siguiente diagrama:

# Ahora, analicemos esta variable por escolaridad

sal_bach<- datos1[escolaridad==0,]$salario

boxplot(sal_bach,horizontal = T)

boxplot(salario~escolaridad,col="yellow") 

# Otra manera
par(mfrow=c(1,1))   # ventana original
library(paletteer)
library(ggthemes)

boxplot(salario[escolaridad=="0"],salario[escolaridad=="1"],
        salario[escolaridad=="2"],salario[escolaridad=="3"],
        col=paletteer_c("ggthemes::Sunset-Sunrise Diverging", 4),names=c("bachillerato","licenciatura sin título","licenciatura con título","posgrado"))

# Otra manera

library(ggplot2)
ggplot(datos1, 
       aes(x = as.factor(escolaridad), 
           y = salario)) +
  geom_boxplot() +
 labs(title = "Salary distribution by rank")

El salario no es el mismo para todos los niveles de escolaridad. Las medianas de los niveles Lst, Lct y posgrado, se parecen. ¿Qué tal si hacemos una prueba de hipótesis para ver esto?

# Prueba NO PARAMETRICA para comparacion de medianas

wilcox.test(salario[escolaridad=="1"],
        salario[escolaridad=="2"], alternative= "two.sided", mu=0, paired= F, conf.level= 0.97)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  salario[escolaridad == "1"] and salario[escolaridad == "2"]
## W = 2712, p-value = 0.9
## alternative hypothesis: true location shift is not equal to 0

De acuerdo con el valor-p se concluye que no hay evidencia estadística suficiente para pensar que la diferencia de medianas no es igual a cero, en otras palabras, la mediana de los salarios para los emplead@s con licenciatura sin título es igual a la correspondiente con licenciatura con título (¡ganas lo mismo con el “papel” que sin él!)

Y si comparamos a los que tienen posgrado:

# Prueba NO PARAMETRICA para comparacion de medianas

wilcox.test(salario[escolaridad=="1"],
        salario[escolaridad=="3"], alternative= "two.sided", mu=0, paired= F, conf.level= 0.97)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  salario[escolaridad == "1"] and salario[escolaridad == "3"]
## W = 92, p-value = 0.6
## alternative hypothesis: true location shift is not equal to 0

¡Ufff! A la empresa no le importa que tengamos un posgrado porque me pagará igual que a la persona que sólo tiene el 100% de créditos de la carrera.

Veamos esto:

library(ggplot2)

# Distribucion del salario 

Nivel.educacion<- as.factor(escolaridad)
num.empleado<- seq(1:200)
ggplot(datos1, aes(x = num.empleado, 
                     y = salario, 
                     color=Nivel.educacion)) +
  geom_point() +
  labs(title = "Salario de acuerdo a la escolaridad")

# O bien,

colors <- c("red", "green", "blue", "orange")
plot(num.empleado,salario,col = colors[Nivel.educacion], 
     pch = as.numeric(Nivel.educacion))
legend("bottomright", legend = c("B", "Lst","Lct","P"),
       lwd = 3, col = c("red", "green", "blue", "orange"))

De acuerdo a las gráficas anteriores, no se observan grupos salariales.

Más sobre el salario:

plot(antiguedad, salario, col=colors[Nivel.educacion], pch = 16)
legend("bottomright", legend = c("B", "Lst","Lct","P"),
       lwd = 3, col = c("red", "green", "blue", "orange"))

# calculate mean salary for each rank
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:Hmisc':
## 
##     src, summarize

## The following object is masked from 'package:MASS':
## 
##     select

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

plotdata <- datos1 %>%
  group_by(escolaridad) %>%
  summarize(sal_prom = mean(salario))

# plot mean salaries
ggplot(plotdata, 
       aes(x = escolaridad, 
           y = sal_prom)) +
  geom_bar(stat = "identity")

library(scales)
ggplot(plotdata, 
       aes(x = factor(escolaridad,
                      labels = c("Bachillerato",
                                 "Licenciatura sin título",
                                 "Licenciatura con título","Posgrado")), 
                      y = sal_prom)) +
  geom_bar(stat = "identity", 
           fill = "cornflowerblue") +
  geom_text(aes(label = dollar(sal_prom)), 
            vjust = -0.25) +
  scale_y_continuous(breaks = seq(0, 50000, 10000), 
                     label = dollar) +
  labs(title = "Salario promedio por nivel educativo", 
       subtitle = "Empresa X",
       x = "",
       y = "")

Para estudiar la relación entre una variable que agrupa (cualitativa) y una cuantitativa podemos usar: strip plot

# Graficamos la distribucion de los salarios por nivel educativo  usando strip plots

ggplot(datos1, 
       aes(y = escolaridad, 
           x = salario)) +
  geom_point() + 
  labs(title = "Distribución de los salarios por nivel educativo")

Es difícil comprender lo que sucede, mejor:

library(scales)
ggplot(datos1, 
       aes(y = factor(escolaridad,labels=c("Bachillerato","Licenciatura\n sin título","Licenciatura\ncon título","Posgrado")), 
           x = salario,color = Nivel.educacion)) +
  geom_jitter() + 
  scale_x_continuous(label = dollar) +
  labs(title = "Distribución de los salarios por nivel educativo")+
  labs(y= "Educación", x = "Salario (USD)")+
  theme_minimal() 

1.5.2 Variable cuantitativa vs variable cuantitativa: ASOCIACION ENTRE VARIABLES

round(cor(datos1[c(1,2,4:6,8:9)]),digits = 2) 

##             antiguedad horasextra cursos incapacidad aptitudes salario edad
## antiguedad        1.00       0.07   0.00        0.05      0.07    0.63 0.76
## horasextra        0.07       1.00  -0.02        0.10      0.14    0.06 0.05
## cursos            0.00      -0.02   1.00       -0.13     -0.10    0.01 0.05
## incapacidad       0.05       0.10  -0.13        1.00     -0.03    0.00 0.04
## aptitudes         0.07       0.14  -0.10       -0.03      1.00    0.16 0.09
## salario           0.63       0.06   0.01        0.00      0.16    1.00 0.75
## edad              0.76       0.05   0.05        0.04      0.09    0.75 1.00

Observamos que las relaciones más débiles se dan entre: antiguedad-horas extra; antiguedad-cursos; antiguedad-incapacidad; antiguedad-aptitudes, horasextra-cursos, cursos-salario, cursos-edad, incapacidad-aptitudes, incapacidad-edad, edad-aptitudes.

De hecho, estas parejas de variables tienen coeficientes cercanos a cero. Así que podríamos sospechar que no hay relación lineal entre ellas. Para verificar esto, hagamos una prueba de hipótesis:

cor.test(antiguedad, horasextra, alternative="two.sided", method="pearson")

## 
##  Pearson's product-moment correlation
## 
## data:  antiguedad and horasextra
## t = 1, df = 198, p-value = 0.3
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.065  0.211
## sample estimates:
##    cor 
## 0.0744

cor.test(antiguedad, cursos, alternative="two.sided", method="pearson")

## 
##  Pearson's product-moment correlation
## 
## data:  antiguedad and cursos
## t = -0.04, df = 198, p-value = 1
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.141  0.136
## sample estimates:
##      cor 
## -0.00266

cor.test(antiguedad, incapacidad, alternative="two.sided", method="pearson")

## 
##  Pearson's product-moment correlation
## 
## data:  antiguedad and incapacidad
## t = 0.7, df = 198, p-value = 0.5
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.0895  0.1873
## sample estimates:
##    cor 
## 0.0499

cor.test(antiguedad, aptitudes, alternative="two.sided", method="pearson")

## 
##  Pearson's product-moment correlation
## 
## data:  antiguedad and aptitudes
## t = 1, df = 198, p-value = 0.3
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.0647  0.2112
## sample estimates:
##    cor 
## 0.0747

cor.test(cursos, horasextra, alternative="two.sided", method="pearson")

## 
##  Pearson's product-moment correlation
## 
## data:  cursos and horasextra
## t = -0.3, df = 198, p-value = 0.8
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.160  0.118
## sample estimates:
##     cor 
## -0.0213

Gráficamente…

plot(datos1[c(1,2,4:6,8:9)])   # matriz de graficas de dispersion

1.5.3 Variable cuantitativa vs variable cuantitativa: COMPARACION ENTRE VARIABLES

Comparamos tanto el salario como las aptitudes:

cv_aptitudes<- sd(aptitudes)/mean(aptitudes); cv_aptitudes   

## [1] 0.118

cv_salario<- sd(salario)/mean(salario); cv_salario    

## [1] 0.222

¿Dónde hay más variabilidad en cuanto al salario de acuerdo con la escolaridad?

bachillerato<- datos1[escolaridad=="0",]  
lic.st<- datos1[escolaridad=="1",]  
lic.ct<- datos1[escolaridad=="2",] 
posgrado<- datos1[escolaridad=="3",] 

(cv_bachillerato<- sd(bachillerato$salario)/mean(bachillerato$salario))

## [1] 0.235

(cv_lic.st<- sd(lic.st$salario)/mean(lic.st$salario))

## [1] 0.222